#########################################################################################
# 
# Main Models 
#
#########################################################################################

#########################################################################################
### Democratic Republic of Congo

m.out.congo <- ictreg.joint(Y ~ female + age.z + edu.z + income.z  + hh_size.z +  murder_yes 
                          + terr_2 + terr_3 + terr_4 + terr_5 + terr_6 + activity_prev,  
                          treat="treatment", 
                          outcome="soc_part.2",
                          outcome.reg="logistic",
                          constrained=TRUE,
                          J=3, data=congo) # List Experiment Model
summary(m.out.congo)

m.out.congo.direct <- glm(soc_part.2 ~ direct + female + age.z + edu.z + income.z + hh_size.z +  murder_yes
           + terr_2 + terr_3 + terr_4 + terr_5 + terr_6 + activity_prev, data=congo, family=binomial("logit")) # Direct Item Model

summary(m.out.congo.direct)

# Post-estimation
sim.eff.fun <- function(m.out){ 
  coef.sim <- mvrnorm(10000, m.out$par.outcome,  m.out$vcov[27:40,27:40]) # Informal Bayes
  
  x.0 <- cbind(m.out$x, rep(0, dim(m.out$x)[1]))
  x.1 <- cbind(m.out$x, rep(1, dim(m.out$x)[1]))
  
  pred.0 <- invlogit(x.0%*%t(coef.sim))
  pred.1 <- invlogit(x.1%*%t(coef.sim))
  
  pred.eff <- pred.1 - pred.0
  pred.eff.av <- apply(pred.eff, 2, mean)
  return(pred.eff.av)
}

sim.eff.fun.2 <- function(m.out){ 
  coef.sim <- mvrnorm(10000, m.out$coefficients,  vcov(m.out)) # Informal Bayes
  
  x.0 <- as.matrix(m.out$model)
  x.0[,1] <- 1
  x.0[,2] <- 0
  
  x.1 <- as.matrix(m.out$model)
  x.1[,1] <- 1
  x.1[,2] <- 1
  
  pred.0 <- invlogit(x.0%*%t(coef.sim))
  pred.1 <- invlogit(x.1%*%t(coef.sim))
  
  pred.eff <- pred.1 - pred.0
  pred.eff.av <- apply(pred.eff, 2, mean)
  return(pred.eff.av)
}


civic.congo.ind <- sim.eff.fun(m.out.congo)
civic.congo.dir <- sim.eff.fun.2(m.out.congo.direct)

#########################################################################################
### Liberia

m.out.liberia <- ictreg.joint(Y ~ female + age.z + edu.z + income.z + hh_size.z + cw_kill +  as.factor(county), 
                          treat="treatment", 
                          outcome="outcome_ca",
                          outcome.reg="logistic",
                          constrained=TRUE,
                          J=3, data=liberia) # List Experiment Model
summary(m.out.liberia)


m.out.liberia.direct <- glm(outcome_ca ~ direct.1 + female + age.z + edu.z + income.z + hh_size.z +  cw_kill 
           + county.1 + county.2, data=liberia, family=binomial("logit")) # Direct Item Model

summary(m.out.liberia.direct)

# Post-estimation
sim.eff.fun <- function(m.out){ 
  coef.sim <- mvrnorm(10000, m.out$par.outcome,  m.out$vcov[19:28,19:28]) # Informal Bayes
  
  x.0 <- cbind(m.out$x, rep(0, dim(m.out$x)[1]))
  x.1 <- cbind(m.out$x, rep(1, dim(m.out$x)[1]))
  
  pred.0 <- invlogit(x.0%*%t(coef.sim))
  pred.1 <- invlogit(x.1%*%t(coef.sim))
  
  pred.eff <- pred.1 - pred.0
  pred.eff.av <- apply(pred.eff, 2, mean)
  return(pred.eff.av)
}

sim.eff.fun.2 <- function(m.out){ 
  coef.sim <- mvrnorm(10000, m.out$coefficients,  vcov(m.out)) # Informal Bayes
  
  x.0 <- as.matrix(m.out$model)
  x.0[,1] <- 1
  x.0[,2] <- 0
  
  x.1 <- as.matrix(m.out$model)
  x.1[,1] <- 1
  x.1[,2] <- 1
  
  pred.0 <- invlogit(x.0%*%t(coef.sim))
  pred.1 <- invlogit(x.1%*%t(coef.sim))
  
  pred.eff <- pred.1 - pred.0
  pred.eff.av <- apply(pred.eff, 2, mean)
  return(pred.eff.av)
}

civic.liberia.ind <- sim.eff.fun(m.out.liberia)
civic.liberia.dir <- sim.eff.fun.2(m.out.liberia.direct)

#########################################################################################
### Sri Lanka

m.out.sri <- ictreg.joint(Y ~ female + age.z + edu.z + income.z + hh_size.z +  killed + trauma + as.factor(Province) + prior,
                          treat="treatment", 
                          outcome="soc_part.2",
                          outcome.reg="logistic",
                          constrained=TRUE,
                          J=3, data=sri) # List Experiment Model
summary(m.out.sri)

m.out.sri.direct <- glm(soc_part.2 ~ direct.1 + female + age.z + edu.z +  income.z + hh_size.z +  killed + trauma +
             + prov.2 + prov.3 + prov.4 + prov.5 + prov.6 + prov.7 + prov.8 + prov.9 + prior, data=sri, family=binomial("logit")) # Direct Item Model

summary(m.out.sri.direct)

# Post-estimation
sim.eff.fun <- function(m.out){ 
  coef.sim <- mvrnorm(10000, m.out$par.outcome,  m.out$vcov[35:52,35:52]) # Informal Bayes
  
  x.0 <- cbind(m.out$x, rep(0, dim(m.out$x)[1]))
  x.1 <- cbind(m.out$x, rep(1, dim(m.out$x)[1]))
  
  pred.0 <- invlogit(x.0%*%t(coef.sim))
  pred.1 <- invlogit(x.1%*%t(coef.sim))
  
  pred.eff <- pred.1 - pred.0
  pred.eff.av <- apply(pred.eff, 2, mean)
  return(pred.eff.av)
}

sim.eff.fun.2 <- function(m.out){ 
  coef.sim <- mvrnorm(10000, m.out$coefficients,  vcov(m.out)) # Informal Bayes
  
  x.0 <- as.matrix(m.out$model)
  x.0[,1] <- 1
  x.0[,2] <- 0
  
  x.1 <- as.matrix(m.out$model)
  x.1[,1] <- 1
  x.1[,2] <- 1
  
  pred.0 <- invlogit(x.0%*%t(coef.sim))
  pred.1 <- invlogit(x.1%*%t(coef.sim))
  
  pred.eff <- pred.1 - pred.0
  pred.eff.av <- apply(pred.eff, 2, mean)
  return(pred.eff.av)
}

civic.sri.ind <- sim.eff.fun(m.out.sri)
civic.sri.dir <- sim.eff.fun.2(m.out.sri.direct)

#########################################################################################
